

# ------------------------------------------------------------------------------
# Code for Shamatha Project 7-Year Follow-Up Study of Adaptive Functioning
# Prepared 2024-04-01

# Load workspace S1.RData by double-clicking on it or opening it from within
# R's native console or your preferred R IDE

# This workspace includes 4 data.frames:
# dL = original data in long format   
# dLi = original data + imputed values (long format)
# dW = original data in wide format    
# dWi = original data + imputed values (wide format)

# Demographic and other identifiable information has been removed.

# Long-format variables are:
# ID					participant IDs, re-numbered for anonymization
# Variable				the measures included in the current analyses
# FixedTime				measurement occasion, where
#						1 = pre-training
#						2 =  0 months post-training
#						3 =~ 6 months post-training
#						4 =~ 1.5 years post-training
#						5 =~ 7 years post-training
# YearPostR				time post-training in years
# PostR					flag marking training vs. post-training periods
# Score					scores on the corresponding variables
# ------------------------------------------------------------------------------

# set random seed for replication purposes
set.seed(613) 

# load libraries
library(lme4)
library(psych)


# ------------------------------------------------------------------------------
# Run the Multilevel Models
# ------------------------------------------------------------------------------

# function to check model convergence
chk_cnvgr <- function(fit, elim_sn = FALSE)
    {
	optmsg  = fit@optinfo$conv$lme4$messages
	nc = ifelse(length(grep("failed to converge", optmsg)) > 0, TRUE, FALSE)	
	sn = ifelse(length(grep("isSingular", optmsg)) > 0, TRUE, FALSE)
	if (elim_sn == TRUE)
	    {
	    prblm = ifelse((nc==TRUE || sn==TRUE), TRUE, FALSE)
		}
	else
	    {
		prblm = ifelse((nc==TRUE), TRUE, FALSE)
		}
	list("prblm" = prblm, "msg" = paste("nc=", nc, ";sn=", sn, sep = ''))
	}

# function for running and comparing models
runmodels <- function(dat, sigcut = .050)
    {
	# fit models (REML by default, which we want for param ests)
    M1 <- lmer(Score ~ PostR + YearPostR + (1 + PostR + YearPostR | ID), 
	            data = dat)   				
    M2 <- lmer(Score ~ PostR + YearPostR + (1 + PostR             | ID), 
	            data = dat)    				
    M3 <- lmer(Score ~ PostR + YearPostR + (1         + YearPostR | ID), 
	            data = dat)   
    M4 <- lmer(Score ~ PostR + YearPostR + (1                     | ID), 
	            data = dat)	
	# test model fits (will auto-refit with ML, which is correct here)
    t1 = anova(M1, M2)
    t2 = anova(M1, M3)			
	t1_x2 = round(t1$Chisq[2], 1); t1_pr = round(t1["Pr(>Chisq)"][2,1], 3)
	t2_x2 = round(t2$Chisq[2], 1); t2_pr = round(t2["Pr(>Chisq)"][2,1], 3)  
	# flag non-convergence and (optionally) boundary singularity warnings
    mE1_cnv = chk_cnvgr(M1)
    mE2_cnv = chk_cnvgr(M2)
	mE3_cnv = chk_cnvgr(M3)
    # choose best model (based on fit comparison & convergence)
    msel = NA					
    # were both random slopes sig/convergent, and was full model convergent?
	if ((t1_pr <= sigcut & t2_pr <= sigcut) & 
	    (mE1_cnv["prblm"]==FALSE) & 
		(mE2_cnv["prblm"]==FALSE) & 
		(mE3_cnv["prblm"]==FALSE)) 
	    {
        msel <- M1; best = "M1"
	# if both random slopes were sig, but full model was non-convergent
    } else if ((t1_pr <= sigcut & t2_pr <= sigcut) & (mE1_cnv["prblm"]==TRUE)) {
	    # check which random slope model had better fit
        if((t1_x2 > t2_x2) & (mE2_cnv["prblm"]==FALSE)) {
		    msel <- M2; best = "M2"	
		} else if (mE3_cnv["prblm"]==FALSE) {
			msel <- M3; best = "M3"
		} else {
			msel <- M4; best = "M4"
		}
	# if both slopes were not sig, figure out which was
	} else if ((t1_pr <= sigcut) & (mE1_cnv["prblm"]==FALSE)) {
        msel <- M2; best = "M2"	
	} else if ((t2_pr <= sigcut) & (mE2_cnv["prblm"]==FALSE)){
        msel <- M3; best = "M3"
	# if both random slopes were non-sig or non-convergent/singular
	} else {
        msel <- M4; best = "M4"
	}
	# model selection info
	sel_info = paste("M1: p=",     format(round(t1_pr,3),nsmall=3), ";",  
		              mE1_cnv["msg"], 
	                 "    M2: p=", format(round(t2_pr,3),nsmall=3), ";",  
		              mE2_cnv["msg"],
					 "    M3: ",                                        
		              mE3_cnv["msg"]
                    , sep = ''
					)
	# grab parameter estimates, individual random effects from best
	results = summary(msel); ranefs = ranef(msel)
	list(c(t1_x2, t1_pr), c(t2_x2, t2_pr), results, ranefs, best, sel_info)
    }

# run 'em
vars = unique(dLi$Variable)											
modres = lapply(vars, function(x) 
    { 
	d = dLi[dLi$Variable == x, ]
    runmodels(d)
    })	
names(modres) = vars


# ------------------------------------------------------------------------------
# Table 1   MLM Results
# ------------------------------------------------------------------------------

# Merge the MLM summary results into a single data frame
tabres <- function(x)
    {
    PrePost_X2    = x[[1]][1]
    PrePost_p     = x[[1]][2]
    Followup_X2   = x[[2]][1]
    Followup_p    = x[[2]][2]	
	fixed = x[[3]]$coefficients	
	randm = x[[3]]$varcor[[1]]
	cors  = attributes(randm)
	resi  = attributes(x[[3]]$varcor)[[2]]
	Ran_Intrcp = randm[[1]]
	Ran_PrePos = Ran_PosSlp = Cor_IPrePos = Cor_IPosSlp = Cor_PrePosSlp = NA
	if (x[[5]] == "M1")
	    {	
		Ran_PrePos    = randm[[5]]
		Ran_PosSlp    = randm[[9]]
        Cor_IPrePos   = cors[[4]][2] 
        Cor_IPosSlp   = cors[[4]][3] 
        Cor_PrePosSlp = cors[[4]][6] 
		} else if (x[[5]] == "M2") {
		Ran_PrePos    = randm[[4]]
        Cor_IPrePos   = cors[[4]][2] 	
		} else if (x[[5]] == "M3") {
		Ran_PosSlp    = randm[[4]]
        Cor_IPosSlp   = cors[[4]][2] 		
		}
		format(round(mean(x, na.rm = T),2), nsmall=2)
		
	res = data.frame("Fix_Intrcp"     = paste(format(round(fixed[1,1], 3), 
		                                             nsmall=3),	
                                              " (", format(round(fixed[1,2], 3), 
                                              	           nsmall=3), ")",
					                          sep = ''),
                     "Fix_PrePos"     = paste(format(round(fixed[2,1], 3), 
                     	                             nsmall=3),		
                                              " (", format(round(fixed[2,2], 3), 
                                              	           nsmall=3), ")",
					                          sep = ''),
                     "Fix_PosSlp"     = paste(format(round(fixed[3,1], 3), 
                     	                             nsmall=3),		
                                              " (", format(round(fixed[3,2], 3), 
                                              	           nsmall=3), ")",
					                          sep = ''),
		             "Ran_Intrcp"     = format(round(Ran_Intrcp, 3), nsmall=3),
                     "Ran_PrePost"    = paste("[", format(round(PrePost_X2, 1), 
                     	                                  nsmall=1), "] ",
		                                      format(round(Ran_PrePos, 3), 
		                                      	     nsmall=3),	
											  sep = ''
											  ),
                     "Ran_PosSlp"     = paste("[", format(round(Followup_X2, 1), 
                     	                                  nsmall=1), "] ",
		                                      format(round(Ran_PosSlp, 3), 
		                                      	     nsmall=3),	
											  sep = ''
											  ),											  
		             "Residual"       = format(round(resi, 3), nsmall=3),		 
					 "Select"         = x[[5]],
					 "MdlSel"         = x[[6]]
                     )
    }			   

# format MLM results	
results <- lapply(modres, function(x) 
    {
	# suppress NA alerts, we're only formatting the data
	suppressWarnings(tabres(x))
    })  	
results <- do.call("rbind", results)	
results$Variable = rownames(results)
rownames(results) = c(1:nrow(results))
ivars = c(11,1,5,13,14,12,4,3,7,9,10,2,6,16,8,15)
cnvrg = data.frame(results[ivars,c(10,8,9)])
results = results[ivars,c(10,8,1:7)]

# TABLE 1
results


# ------------------------------------------------------------------------------
# Table S1.3    Follow-up tests of pre vs. 7-year differences
# ------------------------------------------------------------------------------

# convenience function for creating formulae w/diff. covariates
make_frmla <- function(DV, IVs)
    {
	as.formula(paste(DV, paste(IVs, collapse = " + "), sep = " ~ "))
	}

# function to test pre vs. follow-up differences using regression
test_dif <- function(dat, covs=NULL)
    {
	dat$dif = dat$F3 - dat$Pre
	D = format(abs(round(mean(dat$dif, na.rm=T) / sd(dat$dif, na.rm=T),3)), 
		       nsmall=3)
	if (!is.null(covs))
	    {
		dat$Female = as.factor(dat$Female)
		dat$Retreat = as.factor(dat$Retreat)
		frml = make_frmla("dif", covs)
	    fit = summary(lm(frml, 
		                 data = dat, 
						 contrasts = list(Female = "contr.sum"
								          , Retreat = "contr.sum"
								          )
						))
	    coefs = format(round(fit$coefficients, 3), nsmall=3)
		coefs
    } else {	
	    fit = summary(lm(dif ~ 1, data = dat))
	    coefs = format(round(fit$coefficients[1,], 3), nsmall=3)
	    names(coefs)[3:4] = c("t", "p")
	    coefs = c("Est.(S.E.)" = paste(coefs[1], "(", coefs[2], ")", sep = ''),
	    			coefs[3:4],
	    			"D" = D
	    			)
    }
    coefs
	}

# vars in order of testing
ivars = c(11,1,5,13,14,12,4,3,7,9,10,2,6,16,8,15)
vars = unique(dWi$Variable)[ivars]

# test difs, no covariates
dif_rslts = lapply(vars, function(x) test_dif(dWi[dWi$Variable == x, ]))
names(dif_rslts) = vars
dif_rslts = data.frame(do.call("rbind", dif_rslts))
dif_rslts$p = round(p.adjust(dif_rslts$p, method = "fdr"), 3)

# TABLE S1.3
dif_rslts

## tests w/covariates (MS Table 2) required identifiable data & so are not 
## included here


# ------------------------------------------------------------------------------
# Compile individual random fx scores for Exploratory Factor Analysis
# ------------------------------------------------------------------------------

ranefs <- lapply(modres, function(x) 
    {
    TS = x[[4]]$ID				    	
    TS$ID = rownames(TS)
	nvarTS = ncol(TS)
	lstvar = nvarTS-1
	TS = TS[,c(nvarTS, 1:lstvar)]
    names(TS)[2] = "Intercept"
	nmTS = names(TS)[-c(1:2)]
	if(ncol(TS) == 3)
	    {
		if (names(TS)[3] == "PostR") { TS$YearPostR = NA }
		if (names(TS)[3] == "YearPostR")
    	    {
			TS$PostR = NA 
			TS = TS[,c(1,2,4,3)]
			}
		}
	if(ncol(TS) == 2)
	    {
		TS$YearPostR = TS$PostR = NA
		}
	TS
	})	

varnames = names(ranefs)
for (i in c(1:length(varnames)))
    {
	ranefs[[i]]$Variable = varnames[i]
	}
	
ranefs = do.call("rbind", ranefs)	
ranefs = ranefs[,c(1,5,2:4)]
ranefs[,3:5] = round(ranefs[,3:5], 3)
rownames(ranefs) = c(1:nrow(ranefs))

# limit to variables with significant random effects for changes
# and also limit to random effects where warning of a boundary
# singularity (model converged but ranef was tiny) was not evident
rpp = ranefs[ranefs$Variable %in% c("Extraversion"       ,
                                    "Agreeableness"      ,
                                    "Conscientiousness"  ,
                                    "Neuroticism"        ,
                                    "Openness"           ,
                                    "Mindfulness"        ,
                                    "AvoidantAttachment" ,
                                    #"AnxiousAttachment"  ,
                                    "Diff.Emo.Reg"       ,
                                    "EgoResilience"      ,
                                    #"Empathy"            ,
                                    #"Anxiety"            ,
                                    "Depression"         ,
                                    "WellBeing"          ,
                                    "Disp.Pos.Emo"       ,
                                    "SelfCompassion"     
									),
						            c(1,2,4)]
rfu = ranefs[ranefs$Variable %in% c("Extraversion"       ,
                                    "Agreeableness"      ,
                                    "Conscientiousness"  ,
                                    #"Neuroticism"        ,
                                    "Openness"           ,
                                    "Mindfulness"        ,
                                    "AvoidantAttachment" ,
                                    "AnxiousAttachment"  ,
                                    "Diff.Emo.Reg"       ,
                                    "EgoResilience"      ,
                                    #"Empathy"            ,
                                    #"Anxiety"            ,
                                    #"Depression"         ,
                                    "WellBeing"          ,
                                    "Disp.Pos.Emo"       ,
                                    "SelfCompassion"     
									),
							       c(1,2,5)]

# convert from long to wide format
long_to_wide <- function(dat,sfx)
    {
	names(dat)[3] = "Score"	
	dat = reshape(dat, timevar = "Variable", idvar = "ID", v.names = "Score", 
		          direction = "wide")
	rownames(dat) = dat[,1]
	names(dat) = gsub("Score.", "", names(dat))	
	names(dat)[-c(1)] = paste(names(dat)[-c(1)], sfx, sep = ".")
    dat
	}
rpp = long_to_wide(rpp, "prepos")
rfu = long_to_wide(rfu, "followup")
identical(rpp$ID, rfu$ID)
rchg_sig = cbind(rpp, rfu[,-c(1)])


# ------------------------------------------------------------------------------
# Table 3   Exploratory Factor Analyses
# ------------------------------------------------------------------------------

# number of factors (Figure S3.3)
scree(rchg_sig[,-c(1)])

# factor extraction (3 factors)
fa3o = fa(rchg_sig[,-c(1)], nfactors=3, rotate="oblimin", fm="minres") 

# TABLE 3
print(fa3o$loadings, cut = 0.0, digits = 3);  
round(fa3o$Phi, 3)

		

